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Non-Markovian effects in the evolution of open quantum systems have recently attracted widespread interest, 
particularly in the context of assessing the efficiency of energy and charge transfer in nanoscale biomolecular 
networks and quantum technologies. With the aid of many-body simulation methods, we uncover and analyse 
an ultrafast environmental process that causes energy relaxation in the reduced system to depend explicitly on 
the phase relation of the initial state preparation. Remarkably, for particular phases and system parameters, 
the net energy flow is uphill, transiently violating the principle of detailed balance, and implying that energy 
is spontaneously taken up from the environment. A theoretical analysis reveals that non-secular contributions, 
significant only within the environmental correlation time, underlie this effect. This suggests that environmental 
energy harvesting will be observable across a wide range of coupled quantum systems. 


In recent years, new trends have emerged in the field of 
open quantum system dynamics, particularly in the general 
area of ultrafast, non-Markovian dynamics 12 for which the 
primary steps of natural photosynthesis provide a fascinat¬ 
ing testing ground. In this area the observation of surpris¬ 
ingly long-lasting coherence in a broad range of photosyn¬ 
thetic pigment-protein complexes (PPCs) EH2 as well as 
other coherence effects ED ED has generated tremendous inter¬ 
est in understanding whether quantum effects may underpin 
the near-unit efficiency of the processes of exciton transport 
and charge separation (see IfTOl for a recent review). Early 
studies of PPCs, such as the Fenna-Matthews-Olson (FMO) 
complex QM2, have already established that a dynamical 
interplay of coherent and dissipative dynamics optimizes tar¬ 
geted exciton transfer. Subsequent work has developed further 
the theme that dissipative quantum dynamics may promote the 
efficiency of tasks in photosynthetic and other organic light¬ 
harvesting materials, with a particular focus on the complex, 
structured environments often found in supramolecular sys¬ 
tems fl8ll26l . Motivated by the pressing need to simulate 
accurately system-environment dynamics beyond the Marko¬ 
vian regime, a range of advanced open-system techniques 
have been developed lfl8l f2Tll27]429l . Of these, the TEDOPA 
(Time Evolving Density with Orthogonal Polynomial Algo¬ 
rithm) method has emerged as one of the most powerful for the 
characterization of transient dynamics mm with rigorous 
error bounds l30l . Due to its unique ability to track the many- 
body entangled state of both the system and its macroscopic 
environment, TEDOPA has shed new light into the mechan¬ 
ics of non-equilibrium open dynamics in a range of molecular 
PPCs, solid state and abstract dissipative systems US] ED [32). 

Here we study a biologically motivated model system, 
namely, a molecular dimer dominating the lowest energy 
exciton in the FMO aggregate, and use finite-temperature 
TEDOPA to explore the dynamics within the correlation 
(memory) time r c of the environment, a regime in which sim¬ 
ple master equation approaches tend to fail li33l f34l . Most 
strikingly, we find that the rate and even direction of en¬ 
ergy transfer becomes dependent on the phase of the initial 
superposition of exciton states, allowing certain preparations 


to extract energy during the non-equilibrium evolution of the 
environment. With the aid of a time local master equation 
(ME), we pinpoint the physical origin of this effect to non¬ 
secular contributions of the non-Markovian dynamics of the 
excitonic system due to the quantum interference of dephas¬ 
ing (transversal) and relaxation (longitudinal) fluctuations. In 
many situations, these effects will average out over very short 
time scales. However, the dipolar coupling between the highly 
absorbing pigments in a densely packed PCC realises an exci¬ 
tonic landscape whose energy splittings are typically such that 
the secular time of the excitonic system can become compa¬ 
rable to the environmental correlation time. As a result non¬ 
secular terms are no longer negligible. Our analysis not only 
allows us to rationalise the trends seen in the TEDOPA data 
across the parameter space but in addition demonstrates how 
the energy harvesting process can occur in a much broader 
range of open quantum systems. Non-secular contributions 
have been recently shown to increase the degree of non- 
Markovianity [35] and have been studied in the context of sin¬ 
gle driven systems, where an enhanced secular time emerges 
as a result of dressing by the external field (36l . Finally, we 
suggest how this effect could be observed in biomolecular or 
artificial devices. 

The model - Let us begin by describing a model dimer sys¬ 
tem as composed by two pseudo-spin-1/2 particles (pigment 
sites a and b with transition frequency oj u j, (see Fig. 0 > cou¬ 
pled via an exchange interaction of strength J. The system 
Hamiltonian reads f~ts = + ‘t°~z + j[fT a + < r b - + cr fl cr*), 

where cry (j = x,y,z) are the standard Pauli matrices. Each 
site is subject to the action of a dephasing (transversal) en¬ 
vironment that produces phase randomization but leaves site 
populations unaffected. We model this local interaction by 
coupling the system to a continuum of harmonic oscilla- 
tors Hb = E; J 0 dkhfklb-(k)bi(k), where the independent 
bosonic operators obey |/;/, bf^ = 6k,i6ij , via the Hamilto¬ 
nian ffi = E, (cH + /) /2 dkgj(k) {bfk) + b] (k)j. Without 
loss of generality we will assume the same dispersion relation 
h(k) and coupling strength g(k) for each site. 

The coherent electronic coupling J leads to the formation 
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of delocalised eigenstates (excitons). We restrict the excitonic 
manifold to the diagonalization of the single-excitation sector 
of ( H S , spanned by \e a g h > = \e) a ® | g) b , \g a e b > = |g> a ® \e) b , 
on the basis of biomolecules usually being subject to weak 
external illumination and/or by doubly excited states typically 
being strongly suppressed l37l . The exciton energies are then 
£ 1,2 = ±y with w 0 = V47 2 + A 2 , corresponding to excitonic 
eigenvectors 


/ I 1 ) \ ( cosi sin § U \e a g b ) \ 

\ |2) J \-sin § cosf /\ \g a e b ) )’ 


(1) 


where A = o> a -u b (co a > a> b ), and 6 denotes the mixing angle, 
defined through the relation tan 6 = 2J/A with 0 < 6 < n/2. 



Figure 1. Schematic representation of a dimeric model system sub¬ 
ject to local dephasing ~y d . Coherent inter-site electronic coupling 
J yields the formation of excitonic eigenstates (states ]1) and |2)). 
For J = 52 cm~ l and optical frequencies u> a and a> b , the excitonic 
splitting cuq falls in the far infrared. Exciton dynamics is dictated by 
transverse (dephasing) and longitudinal (relaxation, marked with a 
red double arrow), environmental processes, whose relative strength 
is governed by the exciton delocalization over sites. In the ultrafast 
regime, these two types of fluctuation may interact, which leads to 
long lasting non-secular effects. 

Using TEDOPA we study the (exact) time evolution 
of an excitonic superposition state of the form | i/r) = 
-4= (e‘t |1) + |2)), with a controllable phase For concrete¬ 
ness, we choose a parameter regime as defined by the model 
system provided by sites 3 and 4 in the seven-site FMO Hamil¬ 
tonian of C. Tepidum as taken from EH EE yielding a mix¬ 
ing angle 6 » 7t/4. The local environment is characterised 
by the smooth part of the experimentally fitted super-ohmic 
spectral function of Adolphs and Renger (AR) ll6ll for the 
FMO complex. The reorganisation energy is 35 cnT 1 and the 
background modes are assumed to be initially in thermal equi¬ 
librium at temperature T. Typical values for the spectral width 
yield r c ~ 600 fs while the system’s secular time is Wq 1 ~ 200 
fs. The time scale over which non-secular effects manifest 
can therefore be greatly enlarged as compared to isolated pig¬ 
ments. 


Figure [2| a) shows the time evolution of the excitonic pop¬ 
ulations at T = 277 K for four different initial superposition 
states with different phase We observe that while for an ini¬ 
tial phase £ = 0 or zr/2 the system relaxes monotonically to 
the equilibrium state, with energy being continuously trans¬ 
ferred from the system into the environment, the behaviour 
for an initial phase of f = n or 3zr/2 yields differs radically in 
the early time evolution, with the population of the high en¬ 
ergy exciton |1) becoming larger (population inversion) than 
the population of the low energy exciton for ~ 100 fs. These 
results suggest that for a phase chosen in one half of the com¬ 
plex plane, energy flows from the environment into the system 
at the early times, although the subsequent evolution is a re¬ 
laxation towards a unique equilibrium state where only the 
lowest energy exciton is populated. 

In Figure |2|b) we set the initial phase to be f = n, allow¬ 
ing for energy extraction from the environment, and study 
the excitonic dynamics for different values of the tempera¬ 
ture. Interestingly, increasing the temperature enhances the 
effect, inducing a stronger population inversion and extending 
the duration of the anomalous dynamics before monotonic re¬ 
laxation sets in. 

Theory - A formal derivation of a ME that ensures evolution 
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Figure 2. (Exact) Time evolution of excitonic populations as calcu¬ 
lated by TEDOPA for a dimer system with - oij = 135 cm -1 and 
J = 52 cnr 1 for representative £ = 0, tt/2, n , 37r/2 phases at T = 277 
K in (A) and for temperatures 0 - 300K following the initial state 
preparation with £ = n in (B). In both cases the system is subject 
to the super-ohmic dephasing background characterized by the AR 
spectral density| 
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under a dynamical semigroup (commonly referred as a Lind- 
blad ME), applied to a non-degenerated system, decouples the 
evolution of coherences and populations in the eigenbasis of 
the Hamiltonian If39l . Therefore it will not capture any initial 
phase effect on the exciton population dynamics. Moreover, in 
a Lindblad ME (or a secular Redfield ME) for a two level sys¬ 
tem, the excitonic populations obey a Pauli ME ]40ll . implying 
monotonicity in the evolution at all times, thence forbidding a 
population inversion as observed in Figs. [2|a) and (b). Both 
statements imply that to gain physical insight about the pop¬ 
ulation inversion witnessed by TEDOPA we need to consider 
a more general microscopic approach. Here, we work in the 
Born approximation through the Redfield equation RT1 [42j 
to analyse bath-mediated couplings between populations and 
coherences. In the exciton basis, we obtain a time-local ME 
leading to a coupled set of differential equations for the evo¬ 
lution of excitonic populations and coherences of the form 

P11 = - r rei (0 sitrdpn + r„ (0 sin 2 Gp22 

- ]^sin( 2 e)\r x ns (i)(p n +p 2 i) - *r-L(0(pi2 — P 21 )] 

P 12 = - i (E\ - E 2 )pi 2 ~ 2F d (0 cos 2 6p\ 2 + 2 Y^ r ( t ) sin 2 8p 21 

1 2 

- -sm 2 6(r re i(t)p l2 + r ex (t)p n ) 

- sin(28) \r x ns (0 (pi 1 - P 22 ) + *rL (f) (pi! -P 22 )] 

(2) 

with P 22 = —P 11 , P 21 = P| 2 and 8 as in eq. (fill. Details of 
the derivation as well as explicit expressions for the coeffi¬ 
cients are given in the appendix. (Small) Lamb-shifts have 
been omitted. This system of equations contains the stan¬ 
dard terms found in its secular approximation, namely, pop¬ 
ulation relaxation, T„,/, and thermal excitation, T ex , which are 

_ w 0 

related by a detailed balance condition, T„ (t) = e '«T rd (t), 
and a pure dephasing term with rate (T d ). We have re¬ 
tained the time-dependence of all rates, as we are interested 
in the early time dynamics. In addition to these contribu¬ 
tions, the full Redfield ME contains additional terms that cou¬ 
ple coherences to populations, and coherences to their com¬ 
plex conjugates. We denote these as non-secular (E^ v ) and 
counter-rotating ll33l (or rapidly varying [02]) (r* r ), respec¬ 
tively. By counter-rotating we understand those E’s having a 
time-dependence which oscillates at twice the frequency u>o 
of the excitonic system (these terms would average to zero in 
a coarse-grained or rotating-wave approximation) ll33l . Non¬ 
secular terms have a time-dependence the contains both rotat¬ 
ing (slow) and counter-rotating (fast) components. 

For the initial state \i//) = -2= ( e^ |1) + |2)), the requirement 
Pit (0) > 0 and a Taylor expansion of the decay rates around 
t = 0; T (t) = T' (0) t + O (0) t 2 leads to the condition 

- y/2sin(28) sin(£ + t) > sir? 8 ^' el ^ (l - e~ (3) 

4 1 ns Iw 

for population inversion to occur. In the absence of non¬ 
secular terms, pn (0) < 0 and population inversion cannot 
occur. The positivity of the right hand side of eq. [3] imposes 


an initial phase £ e [3^/4, ln)\\ for population inversion to 
occur, thereby explaining the results of the TEDOPA simula¬ 
tions. Defining the rate of change of energy in the excitonic 
subsystem as A E(t) = wopn(f). we observe that the popula¬ 
tion inversion requires an initial positive AE(0); consequently 
there is a net increase of the energy in the system that is being 
provided by the non-equilibrium environment. Due to the de¬ 
tailed balance condition on the RHS of Eq. [3] the increase of 
the maximum population inversion with temperature, seen in 
Fig. [2[b), is neatly rationalised by the ME analysis. Finally, 
we see from Eq. [3] that the counter-rotating terms do not play 
a role in generating the population inversion, and we will not 
discuss these terms any further. 

The quantitative agreement of the analysis above with the 
main trends seen in the TEDOPA results suggests that a micro¬ 
scopic understanding can be gained by analyzing the structure 
of the non-secular terms in Eq. ([2]). In the excitonic basis, all 
dynamics arise from the environment interactions, described 
through 'TTj = (geos 6cr z - g sin 8cr x ) ® (bf, + b?), where bk 
denote environmental operators. 

The excitonic coupling to the environment therefore in¬ 
cludes a longitudinal (relaxation) component in addition to the 
transversal (dephasing) noise that lead to interference effects 
and the dynamical coupling of populations and coherences. 
Non-secular terms in the time evolution can no longer be ne¬ 
glected, as the resulting quantum interference will be crucial 
to understand the dynamics in the ultrafast time scale. 

Heuristically, the physical origin of the phase-dependence 
of excitonic transport can be understood by considering 
the transition amplitudes between the initial state \\jj) = 
4= |1) + |2)) and the eigenstates 1 1) , |2). A direct cal¬ 

culation yields (1| \H t \i//) = -A |^sin(0) - e 1 ^ cos(0)| Xf bli and 

(2[ |Hi \\p) = -4= [sin(0) + e'f cos(P)J Xf bl j, where the environ¬ 
ment matrix element is between initial and final bath con¬ 
figurations. These amplitudes can be understood simply as 
arising from two interfering pathways; the amplitude for a 
flip |1) <-» |2) being proportional to sin(0) and the amplitude 
for the populations to remain unchanged, |i) —> |i), which is 
proportional to cos (8)e‘? The key observation is that £. con¬ 
trols whether the interference is constructive or destructive 
for a given transition (conservation of probability ensures that 
the other transition is suppressed or enhanced, accordingly). 
Moreover, by virtue of the sin(28 ) proportionality of the non¬ 
secular terms the mixing angle 8 = n/A maximises the pop¬ 
ulation inversion, corresponding to a maximum interference 
between the longitudinal and transversal components of the 
environment while completely delocalised (8 = n/2) or lo¬ 
calised (8 = 0) states lead to no inversion in the population 
evolution. The interference between population preserving 
and inverting pathways is described in the ME approach by 
the non-secular terms. Those terms are characterised for not 
conserving the energy in the system. The physical framework 
is thus that of an electronic transition in the system being com¬ 
pensated for by the creation/anihilation of a virtual phonon in 
the environment. This is a process that drives the environment 
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out of equilibrium, and consequently must vanish in the longer 
time-scales, thereby setting the transient nature of the popula¬ 
tion inversion and ensuring always the relaxation towards a 
thermal state. 

Our ME analysis demonstrates that the phase-dependent 
transport uncovered by the TEDOPA numerics is not a pe¬ 
culiarity of a specific set of parameters but a rather general 
feature of composite open quantum systems, with the mag¬ 
nitude of the population inversion being sensitive to the en¬ 
vironment spectral function and the mixing angle as deter¬ 
mined by the coherent coupling between sites. Indeed, we 
can reproduce the qualitative features of the effect observed 
by TEDOPA employing the Redfield ME. These results, pre¬ 
sented in the S.I. also show that any violation of positivity 
(a known problem of equations beyond the Born-Markov ap¬ 
proximation EH S3) of the reduced density matrix is sev¬ 
eral orders of magnitude smaller in magnitude than the ob¬ 
served population inversion and this together with the fact 
that TEDOPA generates, by construction, a manifestly posi¬ 
tive and numerically exact many-body density matrix, indicate 
that this effect is real and observable in principle. 

Returning to the example of PPCs, most realistic spec¬ 
tral functions contain, in addition to the smooth background, 
a number of sharp features, usually associated with under¬ 
damped intramolecular vibrational modes Ifl8ll22ll38l l44, : 45'l. 
We tackle this more complicated problem with TEDOPA, and 
examine the influence of two such modes (one of them res¬ 
onant with the excitonic gap) in the phase-dependent popu¬ 
lation inversion. Fig. [3] shows that the maximum inversion is 
significantly enhanced (up to 20 % for T = 277 K and £ = n as 
compared to Fig. [2ja]) by the presence of discrete modes. The 
situation is easily understood via an exact quantum mechani¬ 
cal analysis of the system and the environment in the limit t 
—> 0. With factorised initial conditions, the second derivative 
of the high-energy exciton population at time zero is given by: 


Pn (0) = - sin(20)!R [p ei e 2 m S 0 (oj)(2n(aj) + 1 )dto 

Jo (4) 

+ ^ki ? (2n(w,) + 1 ), 

i 

where n(a>) is the Planck distribution. So (k) is the continuum 
part of the spectral function S (to) and i represent an arbitrary 
finite number of intramolecular modes that couple to the sys¬ 
tem with amplitude gj. As it can also be shown that p \i (0) = 0 
(see SI), Eq. [4]shows that discrete modes will always increase 
the initial rate of energy flow, explaining the enhancement of 
the population inversion seen in Fig. [3] Interestingly, the anal¬ 
ysis shows that the initial enhancement due to discrete modes 
does not depend on the frequency of the modes, i.e. exact 
resonance is not required (frequency only appears through the 
temperature dependence). 

Conclusions - By combining powerful numerical tech¬ 
niques with theoretical model analysis, we have demonstrated 
that electronic coherence in a dimeric system can not only 
quantitatively influence the flow of energy transfer over the 
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Figure 3. TEDOPA numerical results for the time evolution of ex¬ 
citon eigenstate populations on a dimer with to a - to b = 135 cnr 1 
and J = 52 cnr' in the superohmic background charaterised by 
the AR spectral density with two discrete vibrational modes of fre¬ 
quencies 180 cnr 1 and 37 cnr 1 at T = 277K for representative 
£ = 0, n/ 2 , 7r, 'in/2 phases. 


correlation time of the environment, it may even revert the 
direction of the flow and permit transient energy harvesting 
from the surroundings. By virtue of the Redfield equation 
analysis, we were able to predict the electronic and environ¬ 
mental properties that maximise this effect, and also showed 
that - in principle - such effects could be found in a wide 
range of natural dimeric systems. As just one possible ex¬ 
ample, the experimental preparation of different excited state 
coherences might be achievable in molecular systems via po¬ 
larisation control or quantum control techniques in 2D Fourier 
Transform spectroscopies. These methods have already been 
shown to be capable of following dynamics within the typical 
bath correlation times we have consider here (w 50 - 100/s) 
ECHOED HZ)- Finally, we remark that the fundamental pro¬ 
cesses we describe strongly motivate the design of thermal 
energy harvesting in quantum devices. As it is known from 
classical systems, these transient effects need to rectified in or¬ 
der to be used, which presents a rather interesting theoretical 
problem in the context of multi-component, non-equilibrium 
quantum dynamics. 
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Dimer parameters and Spectral function for TEDOPA 
simulations 


TEDOPA simulations were carried out for a dimer system 
with local disorder o> a - a>b = 135 cur 1 and coherent cou¬ 
pling J - 52 cm 1 , yielding a mixing angle of ft - tt/ 4. These 
parameter choice is motivated by an actual natural system as 
they correspond to the two lowest sites energies of the FMO 
photosynthetic complex and their electrostatic coupling ED- 
This model system has been previously considered in the dis¬ 
cussion of long lived excitonic coherence If38l . 

The spectral function used to characterise the system- 
environment interaction in TEDOPA simulations is given by 


S{co) = 


lOOOwV 


■V^T + 4.3u> s e 


9! (lOOOwj + 4.3wf) 

2 

+ ^ g.ojfSiio - oj,). 


(5) 


i= 1 


The parameters for the continuum part are A = 35cm 1 , 
c>\ = 0.57cm 1 and ojo = 1 .9cm 1 as taken from the experi¬ 
mental fitting provided in DU. This smooth background dif¬ 
fers from the standard Drude-Lorentz function in being super- 
Ohmic, with the high frequency part of the spectrum respon¬ 
sible for fast thermalization, as required for efficient energy 
transfer. 

The considered discrete modes have frequencies <x>\ = 
180cm _1 and cj2 = 37 cm ' 1 and couplings g\ = 0.12, g2 = 
0.22, as reported in H1 611481 . 


Master equation derivation 


Our theoretical model is provided by an excitonically cou¬ 
pled dimer (sites a and b) with Hamiltonian 


'Hs = y < + y <?crf + j(cr a + o- b _ + crVf) . 


OJb 

~2 


( 6 ) 


Each site is coupled to a vibrational environment modelled 
as an infinite bath of harmonic oscillators, whose Hamiltonian 
is 'Hb = E; J^° dkgi ( k) bj (k ) b, ( k ). The interaction Hamilto¬ 
nian describing the electron-phonon coupling is considered to 
be of the form 

<Hi=Y j h +/ )£° dkh ‘ (k) i bi(k) + b i (k) ) ■ (7) 


We are interested in the time evolution of the diagonal states 
of 9As, known as exciton states. Moreover, we will ignore ex¬ 
citonic states that correspond to either no excitation on each 
site (Igagb)), or to both sites excited (\e a eb)), where the for¬ 
mer are irrelevant to our discussion and the later are typically 
strongly suppressed in the systems of interest (37) . Hence, 


working on the one excitation sector, spanned by the (local¬ 
ized) states \e a gb) and \g a eb), the excitonic states can be ex¬ 
pressed a coherent superposition of the form 


/ |!> \ / cos § sinf W | e a g b ) \ 

\ |2> / \ -sinf cos f J\ | g a e b ) /’ 


( 8 ) 


with energies £| 2 = ±f V47 2 + A 2 = ±^, where A = 
- (x>b ( u> a > u>b) and tan 6 = 2J/A is the so called mix¬ 
ing angle, a measure of the degree of delocalization resulting 
from the coherent coupling J. The excitonic splitting is then 
given by cuq. At this point it is convenient to redefine the en¬ 
vironmental bosonic operators as b p m = ~^(b a ± bb ), which 
when we impose the restriction of one excitation sector, will 
remove the contribution from the b p operators, thus leaving 
just one set of harmonic oscillators. The complete Hamilto¬ 
nian for the effective global system (excitonic + vibrational 
degrees of freedom) is then Tis + 'kin + 77/, where 


7<s = ^w 0 cr, 

<H B = f dkh(k)b\k)b(k) 


'Hi = ( cos6cr z + sin6cr x ) 


dkg(k) (b(k) + b f (k)) . 


(9) 


For simplicity we have assumed that both sites are subject to 
the same type of environment. 

In order to obtain the Master Equation for the evolution of 
the excitonic system, we move to an interaction picture with 
respect to the Hamiltonian 77 = 77s + 77 B . in which. 


A 1 ( t ) = exp (i'Ht)A exp (-177 1) . (10) 


From now on we omit the superscript “I” on the interaction 
picture operators. Assuming a factorized initial state of the 
form p (f) = ps (t) ® Ph and working in the Born-Markov ap¬ 
proximation, the von-Neumann equation for the evolution of 
the density matrix reads 

p(0 = - f dt 1 Tr B {['HHt),['H I (t 1 ),p s (t)®p B ]]} (11) 

Jo 

which is a time local Master Equation, also known as the Red- 
field equation (33) - The coupling Hamiltonian ///(f) in the 
interaction representation is given by 

77/ = [cov0cr. - sinO(coscjQtcr x - sinojoto'y)^ 

r 00 ( 12 ) 

® J dk(b(k)e~ ,u, ‘ + ^^6'“'), 

where we have redefined the bath operators as the product 
g b(k) and assumed the coupling g to be /-independent. The 
presence of crossed terms involving distinct Pauli operators 
will lead to non-secular contributions in the resulting exci¬ 
tonic master equation. To perform the derivation, the follow- 
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ing identities prove useful 

, r = ( 8 jz 8 ix ~ idjy \ 

' U js + tijy S M ) 

[ft";; <Tj\ = 2i£iJkO"k (13) 

{(Ti,(Tj} = 26ijI 

cr x = cr_ + cr + , cr v = i (<x_ - <r + ) . 

Disregarding Lamb shifts, which are unimportant in our 
context, the time-local Master Equation for the excitonic dy¬ 
namics can be written as follows 


j f P (0 = - Id (t) cos 2 6(p (; t ) - cr z p (t) o%) + sin 2 0 (r+ r (t) o-+p(t) cr+ + T nr (t) cr-p (t) cr_) 

+ r re/ (t) sin 2 6\cr-p(t) cr■+ - {cr + o-_,p(f)}j + T ex (t) sin 2 0 |cr + p (t) i {cr_cr + ,p (f)}j (14) 

- \sin(20) (r^ s (t) (cr z p {t) cr x + cr x p (t) cr z ) - (t) ((T-p ( t ) cr y + o- y p ( t ) cr,)) , 


where the expressions for the different rates are given by 


r d (t) = 2 f dt\Q(t - t \), 

Jo 

C(d = 2 f c/f ie i(w » (,+f ' M e(f-fi), 

Jo 


r (0 = 2 


f 


-«(o»o(?+0)) 


e(f-D, 


E re / (0 


dt\cos (a>o (t - ti)) Q(t - t\), 


r ex (t) = 4 f dt\cos (u>o (t - t\)) e~^Q (t — ti), 
Jo 




- 2 f 

•■f 


aio(t-ti) ojo(t + ti) 

dticos --- cos --- Q (t - 1 1 ), 


T -v r, I ■ ojo(t-ti) . a>o (t + h) 

I Jv (t) = 2 | dt\ sin --- sin --- Q(t-h). 


(15) 


The function Q(t - t\) gathers the effect of the environment 
as 

XOO / Q \ 

dot S (u>) coth^—jcos (u>(t - ti)), (16) 

with the spectral density S(oj) describing the weighted dis¬ 
tribution in frequency space of the environment of harmonic 
oscillators. 

The dynamics resulting from the considered ME approach 
is illustrated in Fig. 4, where we show the time evolution of 
the excitonic populations and the rate of change of energy, as 
well as a positivity test of the resulting (non-secular) density 


matrix. In these results we introduced heuristically a decay 
for the non-rotational and non-secular rates via the time de¬ 
pendent exponential exp{-u>ot}, while for simplicity leaving 
the rest of the rates equal and constant. 


Exact derivation with discrete modes 


Let us analyse the effect of the presence of intramolecular 
vibrational modes by performing an exact quantum mechan¬ 
ical calculation, which is possible at t — 0 when assuming a 
factorized initial condition. Using the excitonic basis for the 
dimeric system, we can write 


E n \n) (n\ + ^ U) k c] k c uk + nm \n){m\. (17) 

i= 1,2 i,k n,m 


Here we have assumed a finite number k of independent bath 
modes coupled to each site of the dimer. In the delocalised 
exciton basis the coupling system-environment gets redefined 
through the new couplings Q nm = 2 i,kgi,kC‘ n C l m (Ci yk + c] k ), 
where the coefficients C can be taken as real, so Q nm = Q mn . 
These coefficients can be re-expressed in terms of the mixing 
angle of the dimer system to obtain explicit expressions for 
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0 0.005 0.01 0.015 0.02 0.025 0.03 

Time 27iQ/to 0 



Figure 4. Master Equation simulation (a) Excitonic eigenstate popu¬ 
lation obtained by numerically solving the ME Eq. [TTJfor represen¬ 
tative £ = 0. tt/2, 7i, 3tt/ 2 at T = 277K. Note that the crossing is qual¬ 
itatively reproduced even in this simple situation, (b) Rate of change 
of energy A£(f) = pn(t) for representative £ = 0,7 t/2,7t, 3tt/ 2 at T 
= 277K, positive the differences when the crossing is present, and 
always asymptotically going to zero as required by returning to equi¬ 
librium. (c) Positivity test for representative £ = 0, tt/ 2, 7r, 3zr/2 at T 
= 277K. Although it is not ensured, note that the scale of violations 
would be orders of magnitude smaller than the excitonic difference 
in the crossing on (a). 


the Q couplings 


G 11 = 

k V 

^)(cu + c}J 

+Y gksinl {[ 

k 

i )( C 2 ,k + c{k) 

Q 22 = X gksin 2 1 ; 

k V ' 

!)(cu + cy 

k ' 

2 ) ( C 27- + c 2 ,k) 

Gl2 = ~2 Y gkSin 

(8) (cijc + c\ k/ 


+ \ Y\ gksin ( 6 ) (c 2 ,k + c\ k ) . 

k 

Introducing the mode displacements Xk = c> + c[ allows us to 
write 


Gn _ G 22 = cosO'^gk (A'l.r- - X 2j k) 

sind v 1 / (19) 

Qll - ~-r- / , gk (X\,k ~X 2 ,k) ■ 

k 

Using the Heisenberg equations of motion 0(t ) = 

—i [*H, 0{f)\ to obtain the equations of motion for the exci¬ 
tonic population and coherence operators we obtain 


Pit (0 = ~iQi 2 (t)(pi 2 ( t ) — P21 (0) 

P 12 (0 = -i (£12 + 611 (0 - G 22 (0)Pi2 ( f ) (20) 

~ iQl 2 (t ) (P22 (0 “Pit (0) . 


Following the analysis performed on the main text, we want 
to know the influence that a phase of a superposition of ex- 
citons has on the direction of the energy flow. The rate of 
change of energy is A(f) = 0 JqP\ \(/), so the sign of pw(t) at 
t - 0 tells in which direction the energy flows after a sharp 
excitation. Thus, taking the second derivative of the high en¬ 
ergy exciton population evolution equation and inserting the 
expression above for the rate of change of coherences we ob¬ 
tain 


P 11 (0 = G 12 (0) (£12 + G 11 (0) - Q 22 (0)) 91 ip 12 (0)). (21) 

Assuming separability into electronic and thermal environ¬ 
ment density matrices, the expectation of the second deriva¬ 
tive of the population with respect to the latter is 

Pi, (0 - 2<G, 2 (0) (Gn (0) - Q 22 (0))pK {p 12 (0)) (22) 

where the thermal expectation value of the environment oper¬ 
ators is 

<Gi 2 (0) (Gn (0) - G 22 m P = -\sinQ® Y + 1). 

k 

(23) 

with n(u>) being the thermal Planck distribution. (Q,j (f)) van¬ 
ishes for a thermal state (explaining why the first derivative 
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W(2«(w,)+ 1) = SM(2«( W )+1)^. (25) 

/ Jo 

And if we now choose the spectral density to contain dis¬ 
crete vibrational modes as S(ta) = Sq(io) + - w,), 

with S o (oj) a smooth function representing the background 
fluctuations and i discrete modes of frequency o), and coupling 
strength g,, we obtain 

/~*oo 

Y j gl(2n(aj k )+]) = S(oj)(2n(oj) + \)doj = 

i Jo 

(26) 

S 0 ((o)(2n(oj) + 1 )doj + ^gf(2n(w,) + 1 ). 


Figure 5. TEDOPA numerical results for the time evolution of exci- 
ton eigenstate populations on a dimer with ot a - at/, = 135 cm~ l and 
J = 52 cm -1 in the presence of two vibrational modes of frequencies 
180 cm -1 and 37 cm~ l for temperatures 0-300 K following the initial 
state with £ = n in the AR background. 


of the populations vanishes as t — 
those two last equations to obtain 


0). We can now combine 


d~p it 
dt 2 


=o = sin (20) % {p n (0)) £ g 2 k (2n(to k ) + 1). (24) 


Introducing now the definition for the spectral function of the 
environment as S (at) = 2,- g 2 6(co - at/) allows us to write: 


Putting together this result and the equation for the second 
derivative of the exciton population Eq- [24] we observe that 
the presence of discrete modes will always enhance the initial 
rate of energy flow, whichever its direction. 

We recover here the result that the phase sensitivity is intro¬ 
duced in the system via a correlation of transversal and lon¬ 
gitudinal components of the environment ((Qn (0) (Qn (0) - 
Qii (0))) in Eq [24] in complete agreement with the ME ap¬ 
proach presented above, where it is also obtained that the op¬ 
timal mixing angle 6 = n/A maximises the effect (overlap) of 
these cross correlations. In addition, the quantum mechanical 
analysis, although valid only at t = 0 (afterwards no longer 
can there be factorisation), shows that the populations grow 
quadratically at early times, as obtained with TEDOPA. 











